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Abstract 

Circadian oscillation provides selection advantages through synchronization to the daylight cy- 
cle. However, a reliable clock must be designed through two conflicting properties: entrainability 
to properly respond to external stimuli such as sunlight, and regularity to oscillate with a precise 
period. These two aspects do not easily coexist because better entrainability favors higher sensitiv- 
ity, which may sacrifice the regularity. To investigate conditions for satisfying the two properties, 
we analytically calculated the optimal phase-response curve with a variational method. Our result 
indicates an existence of a dead zone, i.e., a time during which external stimuli neither advance 
nor delay the clock. This result is independent of model details and a dead zone appears only 
when the input stimuli obey the time course of actual insolation. Our calculation demonstrates 
that every circadian clock with a dead zone is optimally adapted to the daylight cycle. Our result 
also explains the lack of a dead zone in oscillators of mammalian somatic cells, and justifies more 
effective entrainment by dawn and dusk than by on/off lighting protocols in animals including 
humans. 
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I. INTRODUCTION 



Circadian oscillators are prevalent in organisms from bacteria to humans and serve to 
synchronize bodies with the environmental 24 h cycle [l, 2|. Environmental selection pres- 
sures, such as that of harmful ultraviolet (UV) radiation, is assumed to be its evolutionary 
cause , and the selection advantage of having the ability to synchronize has been 
experimentally demonstrated by using photosynthetic species such as cyanobacteria and 
Arabidopsis [5ru|- Interestingly, the molecular implementation of oscillation is species spe- 
cific [8]: only three clock proteins can drive oscillation in cyanobacteria |9J], whereas in 
mammals, there is interplay of multiple genetic feedback loops lOj] . Still, all of them satisfy 
two requirements to achieve the reliable synchronization to the environmental cycle: en- 
trainability to properly respond to external stimuli such as sunlight and regularity to 
oscillate with a precise period. A major source of interference with regularity is discreteness 



of molecular species, i.e., molecular noise 



llNl5|. Many studies have analyzed the resistance 
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3. 



Regarding entrainability, cir- 



mechanisms of circadian oscillators against the noise 
cadian clocks synchronize their internal time with the environmental cycle via sunlight, and 
its effect depends on the wavelength or fluence, as well as on the phase of the stimulation. 
However, entrainability and regularity are conflicting factors, because circadian clocks with 
better entrainability are sensitive not only to the light stimuli, but also to the molecular 
noise which interferes with regularity. 

Since both regularity and entrainability are important adaptive values, we would expect 
actual circadian oscillators to optimally satisfy these two factors. Here we investigate the 
optimal phase-response curve (PRC), which is both entrainable and regular, in the phase 



oscillator model 



191 ] by using the Euler-Lagrange variational method. Our main finding 



is the inherent existence of a dead zone in the PRC: optimality is achieved only when the 
PRCs have a time period during which light stimuli neither advance nor delay the clock 
(Fig. (IK). In other words, a PRC with a dead zone (Fig. Q]A.) is better adapted than those 
without a dead zone (Fig. [I]B). We also tested this with two types of input stimuli: an 
insolation-type input that simulated the time course of solar radiation intensity (cf. Eq. [T2 
and Fig. [2A) and a simple sinusoidal input (sine curve). Surprisingly, the dead zone in the 
optimal PRC only emerges for the insolation-type input, not for the sinusoidal input. Along 
with a fact that many experimental studies reported the existence of the dead zone , our 
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results indicate that circadian oscillators in various species have adapted to insolation to 
improve synchronization. 



II. THEORETICAL BACKGROUND 



Phase oscillator model 



Circadian oscillators basically comprise interaction between mRNAs and proteins, whose 
dynamics can be modeled by differential equations. A circadian oscillator of A-molecular 
species can be represented by 

d ^ = F i {x) (z = l,2,.-.,iV), (1) 
where the A- dimensional vector x = (xi, X2, • • • , x^) denotes the concentration of molecular 



species. By using the phase- reduction approach [19|, |20j, the A-dimensional equations of 
unperturbed limit-cycle oscillators can be represented by = Q = 2tt/T, where Q and T 
are the frequency and the period of the oscillation (0 e [0, 27r)), respectively. The effect of 
noise on genetic oscillators has been a subject of considerable interest, and noise-resistant 



mechanisms have been extensively studied 2lH23|. In general, the dynamics of the i-th 



molecular concentration in a circadian oscillator subject to molecular noise is described by 
the following Langevin equation: 

(IT ■ 

-^ = F l (x;p) + Q l (x)Ut), (2) 

where Qi(x) is an arbitrary function representing the multiplicative terms of the noise, £j(t) 
is white Gaussian noise with the correlation {£,i(t)£,j(t')) = 25^5(t — £') (a bracket (•) denotes 
expectation), and p is a model parameter. 

Circadian oscillators synchronize to environmental cycles by responding to an input sig- 
nal. We let p in Eq. [2] be stimulated by the input signal: for example, p can be the degra- 
dation rate. (For simplicity, we consider that the input signal affects only one parameter.) 
We use Eq. [2] for calculating regularity and entrainability of circadian oscillators. 

B. Definition of regularity 

Because the circadian oscillator of Eq. [2] is subject to noise, its period varies cycle to cycle 
(see the supporting information (SI) appendix, section 1). We use term regularity for the 



period variance of the oscillation. Let us first consider the case without input signals (i.e., 



p is constant). By using standard stochastic phase reduction [l9j, Eq. [2] can be transformed 
into the following Langevin equation with respect to the phase variable (Stratonovich 
interpretation): 

d ^ = n + jru l (m.(<i>Mt), (3) 

i=l 

where U{<f>) = (C/i(0), • • • , U N (<j>)) is an infinitesimal PRC (iPRC) U{<f>) = V x (j>\ x=Xhc(4>) , 
and we abbreviated Q«(cclc (</>)) as Qi{4>). The N- dimensional vector £Clc(</>) denotes a point 
on the limit-cycle trajectory at phase 0, where LC stands for limit cycle. The value of iPRC 
Ui{4>) is calculated as a solution of an adjoint equation |20| or as the set of eigenvectors 



of a monodromy matrix in the Floquet theory [19( for arbitrary oscillators. In Eq. [3j the 
average period corresponds to the mean first passage time with <p starting from to 2ir, 
and the period variance is the variance of the first passage time. Direct calculation of 
the period variance is difficult; we first calculate the phase variance and then approximate 
the period variance by multiplying the phase variance by T 2 /(47r 2 ), after Kori et al. 24]. 



Introducing a slow variable (p = (p — ilt and applying phase averaging [19], the phase variance 
after one period T can be calculated since the time evolution of ip approximately obeys 
a simple diffusion equation. Considering the ratio between and Vr, the period variance 
is approximated as 

/ T \ 2 T 3 f 27T N 

vt-v^-J =^j Q deY^um'Qm 2 - (4) 

Details of the calculations are given in the SI appendix, section 1. 



C. Definition of entrainability 

The entrainment property is an important characteristic of limit-cycle oscillators and 



attracts attention in systems biology 25M28I] . For instance, a period mismatch in coupled 



oscillators is known to enhance entrainability in genetic oscillators [281 ] . Light stimuli affect 
the rate equations, i.e., the parameter p in Eq. [2] is perturbed as p + dp by the input signal. 
Equation [2] can be viewed as representing the dynamics of a tilted periodic potential (i.e., 
ratchet) subject to noise. Since a synchronizable condition corresponds to the existence 
of stable points in the ratchet-like potential, the entrainability can be discussed without 
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considering the noise. Consequently, in contrast to the calculation of regularity, in the 
evaluation of the entrainability, we consider a case without molecular noise (i.e., Qi(x) = 
in Eq.EJ. 

Let p(ut) be an input signal with angular frequency u. Considering a weak periodic input 
signal dp = xp(wt), where x is the signal strength (x > 0), and applying the phase reduction 
approach to Eq. [2J the time evolution of the phase variable <fi is given by = Q + xZ(<p)p(u)t), 
where p) = Fi(x LC {<j))] p) and 

= (5) 



8=1 



Z(0) is the PRC with respect to the parameter p and corresponds to experimentally observed 
PRCs. In order to distinguish Z(<p) from iPRC Ui(<p), we will refer to Z((f>) as the parametric 
PRC (pPRC) 29j. Note that phase reduction can yield reliable results only when the 
perturbed trajectory is close to the unperturbed limit-cycle trajectory (i.e., x is sufficiently 
small). We calculate both iPRC and pPRC in our analysis. 

We next quantify the extent of synchronization to the periodic input signal. By intro- 
ducing another slow variable ip = <f) — ut, we can again apply the phase- averaging procedure, 
which yields ip ~ AQ + xG (■?/;) with AQ = Q — u and 

e(i>) = ±-J\ez(i> + 6) P (9). (6) 

The oscillator of interest can synchronize to external signals when there is a stable solution 
of ip in = Afi + x@(V'), i.e., x@(V'm)+^ < oj < x < S>(4'm) + ^ with ipM = argmax^, 9 (■?/>) and 
ijj m = argmin^B (?/>). We can define the entrainability S, or the extent of synchronization, 

as 

£ = e(i> M )-Q(i> m ), (7) 

where x& constitutes the width of the Arnold tongue for sufficiently small x (see the SI 
appendix, section 2). 

Intuitively, a circadian oscillator with better entrainability (i.e., larger E) can synchronize 
to an input signal that has a period further from that of the oscillator. The calculation above 
is standard in the phase reduction approach, and further details are available in Ref. [ig| ] . 
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D. Variational method 



We use the variational method to calculate the optimal PRCs, which maximize the en- 
trainability £ subject to constant variance Vt = <j\- While the conventional optimization 
method can only search for parameter values, the variational method can obtain the op_ 



timal function shape (i.e., the PRCs). In the context of neuronal oscillators, a study [30 1 
has used the variational method to calculate the optimal PRCs for stochastic synchrony 
(noise- induced synchronization). 

The variational equation to be optimized is 

C[U] = £[U) - XV T [U], (8) 

where A is the Lagrange multiplier. Note that variational Eq. [8] is similar to Ref. which 
optimizes the input signal for the maximal entrainment under constant power of the input. 
The variational condition 5C[U] = yields the optimal iPRC 

and the pPRC is calculated with Eq. Since ipM and ip m themselves depend on t/j(0), 
they have to satisfy a self-consistent condition, i.e., Eq. [7] is maximal with ipM and ip m . 
Consequently, we maximize the following function: 



with 



\/^(A,S), (10) 







Q t .(6 + 5y V dp 
where A = ipM — ipm an d 5 = ip m . The optimal iPRC can be obtained by first finding the 
maximum solution of ^(A, 8) with respect to A and 5, and then substituting the obtained 
solution ip m = S and ipM = 6 + A into Eq. [9] (the pPRC can be calculated by Eq. [5]). 



E. Input signal to model insolation 

Optimal PRCs depend on input signals, as seen in Eq. [91 The most common synchronizer 
in circadian oscillators is light stimuli, for which the strength is determined by 24 h-periodic 
insolation. The insolation is calculated by I = / O cosi? for < $ < 7T, where $ is the zenith 
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angle, Jo is the maximum insolation, and 1 = when the sun is below the horizon 32] . It 
can be approximated by 

piuit) = ramp(sin(u;t)), (12) 

where ramp(x) is the ramp function defined by ramp(x) = x for x > and ramp(x) = for 
x < 0. We call Eq. [12] the insolation input, and the plot is shown in Fig. |2JA. (the shaded 
region represents night). For general discussions on the daily insolation, readers are referred 



to Ref. 



33| (an eBook). For comparison, we also employ a sinusoidal input, which is common 



in nonlinear sciences: 

p(cjt) = 8m(ut). (13) 

For calculating the optimal PRCs, we use Eqs. [12] and [13] (a squared sinusoidal insolation 
input is considered in the SI appendix, section 6). 



III. RESULTS 



A. Optimal PRC shape does not depend on the period, its variance, or noise 
intensity 

In all specs, light stimuli affect the oscillatory dynamics multiphcatively, i.e., they act on 
the kinetic rates or transcriptional efficiency of the gene regulatory circuits |8|]. We assume 
that the j-th molecular species includes a parameter p as 

Fjix; p) = Fj(x) + px k , (14) 

where Fj(x) represents the terms that do not include p, and x k is the concentration of the 
fc-th molecular species. Here, k e {1,2, ••• ,N} can take any value regardless of j (both 
j 7^ k and j = k are allowed). When using phase reduction, the dynamics of the limit 
cycle are considered on the unperturbed limit-cycle trajectories x^c, an d hence the points 
on the limit cycle can be uniquely determined by the phase <fi. Consequently, under the 
phase reduction, x k is replaced by £Lc,fc(0) lYl Eq- UH where XLc,fc(0) is the fc-th coordinate 
of x LC (i.e., d p Fj((p; p) = XLc,fc(0) in Eq. Here, xlc,A;(0) corresponds to the time course 
of the concentration of the k-th molecular species. Because x^c,k{4>) constitutes a core 
clock component and is generally a smooth 27r-periodic function, we approximate it with a 
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sinusoidal function: 

xlcA^) = 1 - asin(<f) + u) , (15) 
where u is the initial phase and a denotes the amplitude of the oscillation (Fig. [2j3) . To 
ensure £LC,fc(0) > 0, we set < a < 1, and a = recovers the additive case. Since the initial 
phase u does not play any role (u is offset by 5 in Eq. ITT1) when the white Gaussian noise is 
additive (i.e., Qi(x) oc 1), we also set u = 0. The parametric approximation of Eq. [TBI enables 
an almost closed form for the overall calculations (see the SI appendix, section 4). For a 
multiplicative noise term Qi(x), we assume that the white Gaussian noise is additive and 
is present only in the j-th coordinate equation (Qj(<f)) = ^/q, where q is the noise intensity 
and Qi(4>) = for i ^ j), because the multiplicative noise model makes overall calculations 
much more complex. 

Figures EA-C show the landscape of \I/(A,<5) as functions of A and 5, and Fig. |3p-F 
express the optimal iPRCs Uj(<fi) and pPRCs Z((fi) for the insolation input. Interestingly, 
the optimal PRCs do not depend on the model parameters such as the period T, its variance 
erf,, or noise intensity q. These three parameters only act on the magnitude of the PRCs (i.e., 
the vertical scaling of the PRCs). Consequently, we normalized T — 1, a\ — 1, and 5 = 1, 
as shown in Fig. |3j As the optimal PRCs depend on a, \I/(A,£) is plotted for three cases: 
(A) a = 0, (B) a = 0.5, and (C) a = 1.0, where the maximal points (A, 5) yield the optimal 
PRCs using Eq. |9j The maximal parameters A and S are numerically calculated with the 
differential evolution. Figures |3p-F describe the optimal iPRCs (solid line) and pPRCs 
(dashed line) for a = 0, 0.5, and 1.0, respectively. When a = 0, i.e., the input signal is 
additive, ^/(A, S) achieves a maximum for A = 7r and arbitrary 5, yielding sinusoidal PRCs 
as the optimal solution (Fig. EP). Although the input signal p(<f)) is not sinusoidal, the 
optimal PRCs obtained using the variational method become sinusoidal. In other words, 
considering optimality, resonator-type oscillators have an advantage over integrator- type 
oscillators. For a > 0, the input signal p(4>) depends on the concentration of the k-th 
molecular species. From Fig.|33, the optimal parameters for a = 0.5 are (A, 5) = (2.67, 1.81) 
and (3.62,4.47), which are different from A = n (these two sets yield symmetric PRCs with 
respect to the horizontal axis). Figure [3E shows the optimal iPRCs Uj(<p) and pPRCs Z(<f)) 
for a = 0.5. Interestingly, the optimal iPRCs and pPRCs for a = 0.5 have a dead zone 
(region of 1 < < 2 in Fig. |3jC) in which the input signal neither advances nor delays the 
clock. From Eq. [9] and the insolation input of Eq. [121 the optimal PRCs inevitably include 
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a dead zone if the optimal A is not tt. For a = 1.0, there are four sets of parameters (A, 5) 
that give optimal PRCs: (2.30,2.72), (2.30,1.26), (3.98,3.56), and (3.98,5.02) (PRCs with 
these four sets are symmetric each other with respect to the horizontal axis or = 37r/2). 
Consequently, the optimal PRCs shown in Fig.[3F have a dead zone as in the case of a = 0.5. 

B. Dead zone is inherent to insolation input 

From the results discussed above, the optimal PRCs have a dead zone when a > 0. We 
next studied the length of the dead zone as a function of a and improvements in the entrain- 
ability induced by the dead zone for the insolation input, as shown in Fig. H]/Y. Figure H]A. is 
a dual-axis plot of the length of the dead zone (solid line; left axis) and the improvements 
in the entrainability (dashed and dot-dashed lines; right axis) as a function of a. Because 
the dead zone emerges from Eqs. M and [J2] when the optimal parameter is A ^ tt, we can 
naturally define its length as 



where A is the maximum value of \&(A, 5). As seen in Fig. H]A., a dead zone clearly exists 
when a > 0, and the length increases with increasing a for a < 0.8. Even for a = 0.1, when 
the oscillation amplitude of £lc,A:(0) (the concentration of a molecular species modulated 
by the light-sensitive parameter p. cf. Fig. [2B) is very small, we observe a dead zone with a 
length of L = 0.475, which corresponds to about 3 h within 24 h, indicating the universality 
of having a dead zone in order to attain optimality. The improvement in the entrainability 
that is induced by a dead zone is calculated by comparing the entrainability of the optimal 
PRCs with that of typical sinusoidal PRCs. We consider sinusoidal functions for both the 
iPRC Uj{(t>) and pPRC Z{<f>) by setting 



where c is the parameter to be optimized so that entrainability is maximized for each a (see 
the SI appendix, section 5 for the explicit expressions). Equations [T7I and [TBI are scaled so 
that they satisfy the constraints on the period variance (Eq. Hj). We calculated the ratios 



L = |A - tt|, 



(16) 



Uj(4>) oc sin(0 + c), 
Z((f>) oc sin(0 + c), 



(17) 
(18) 



Ru — £/£u, Rz — £/£z 



(19) 
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where 8u and £z represent the entrainabilities for the cases of the sinusoidal iPRC and 
pPRC, respectively, which were calculated with insolation input. For the sinusoidal iPRC 
of Eq. [T71 the entrainability is calculated with pPRC via Eq. [5j Ru and Rz quantify the 
improvement rate of the optimal PRCs over the sinusoidal iPRC (Ru) and pPRC (Rz)- In 
Fig. HA, the dashed and dot-dashed lines show Ru and Rz, respectively, as a function of 
a. Both ratios monotonically increase as a increases, which shows that the optimal PRC 
with a dead zone exhibits better entrainability when the oscillation of XLC,fc(0) has a larger 
amplitude. When the concentration of £lc,Ai(0) is low, the effects of the input signal on the 
circadian oscillators are smaller. This is because pPRC Z(<f>), which quantifies the extent 
of the phase shift due to the stimulation of the parameter, depends on the concentration 
x LC,k(4 ) ) (see Eq. [5]). However, even within the range where £LC,fc(0) has smaller values, 
the iPRC Uj(4>) contributes to an increase in the variance of the period, regardless of the 
concentration. From this, we see that having an iPRC with a smaller magnitude when the 
concentration of £lc,A;(0) is smaller results in a smaller variance, which results in a larger 
entrainability for a constant variance of the period. Figure [5] shows dual-axis plots of the 
iPRC (solid line; left axis), the pPRC (dashed line; left axis), and the concentration xlc, &(</>) 
(dot-dashed line; right axis) for (A) the optimal case, and (B) an alternative case that has a 
dead zone in the range in which the concentration of x^c,k(4 > ) nas a peak (i.e., opposite from 
the optimal case). Since both cases share the same iPRCs Uj(<fi) except for a difference in 
phase, these two iPRCs yield the same variance in the period. However, the corresponding 
pPRCs due to Eq. [5] are different because the concentration xlc, &(</>) of the part that can be 
phase shifted is higher for the optimal case, which results in better entrainability. Figure H]B 
compares the pPRCs Z(<f>) with a = 0.5 for three cases: the optimal pPRC identical to 
the dashed line in Fig. [3] (solid line), the pPRC of the sinusoidal iPRC Uj(<f>) according to 
Eq. [17] (dashed line), and the sinusoidal pPRC Z(<j>) of Eq. [18] (dot-dashed line). These three 
pPRCs Z(4>) have corresponding iPRCs Uj(<p) that yield the same period variance in Eq. 0] 
Although the examples in Fig. [5] qualitatively explain the benefits of a dead zone, for some 
input values, the optimal PRCs may not contain a dead zone for any value of a. This will 
be shown next. 
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C. Sinusoidal input does not result in a dead zone 

Since the optimal PRCs depend on input signals (Eq. [9]), we next consider a typical 
periodic input signal, a sinusoidal function (Eq. [T3]). In this case, ^(A, S) is calculated in 
a closed form (see the SI appendix, section 4), and ^(A, S) yields the maximal value for 
(A, 5) = (tc, rnr) for < a < 1, where n is an integer and when a = 0, 6 can take any value. 
The optimal PRCs due to Eq. [9] do not exhibit a dead zone (Fig. [6]). For a = 0, the optimal 
PRC is sinusoidal and is equivalent to Fig. [3p. For a = 0.5, the optimal PRC is still close to 
a sinusoidal function (Fig. [6A). When increasing a to a = 1.0, the PRC diverges from the 
sinusoidal function and exhibits almost positive values (Fig. EB); in neuroscience research, 
this is called an integrator-type oscillator. 



IV. DISCUSSION 



The existence of a dead zone optimizes both entrainability and regularity. It is obvious 



that having a dead zone optimizes stabi 
any kind of signals, which was shown in 



ity alone, as such null PRCs are never affected by 



34j with respect to daylight fluctuations. However, 



optimality of both entrainability and regularity, which are in a trade-off relationship, con- 
ferred by a dead zone is highly non-trivial. Furthermore, our finding is fairly general since 
it depends very little on the details of the model. The optimal PRCs are determined by A 
and 5 via Eq. O which were obtained as the parameters that maximize \I/(A,5). Although 
there are four parameters in our model (T: period, a: period variance, q: noise intensity, 
and a: oscillation amplitude), A and S depend only on a. Furthermore, as seen in Fig. HJ\, 
a dead zone always exists in an optimal PRC unless a = (additive stimulation). We also 
calculated the optimal PRCs for a power of a sinusoidal function, and obtained a result 
for a > that was qualitatively the same (see the SI appendix, section 6). Along with 
the fact that T, <tt, and q affect only the scaling of the optimal PRCs, when the input 
signal affects the dynamics multiplicatively (i.e., a > 0), the existence of a dead zone always 
provides a synchronization advantage. This is supported by many experimental studies of 
various species, that report the existence of a dead zone in the PRC |l|. Our general result 
suggests that circadian oscillators have fully adapted to insolation to improve synchroniza- 
tion. Indeed, many experimental findings imply that circadian oscillators have adapted to 
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actual insolation 



35[: for various animals, including scorpions, fish, hamsters, monkeys, and 



humans, light-dark (LD) cycles that include a twilight period result in better entrainability 
than do abrupt LD cycles (on-off protocols) 35]. 

The insolation input plays an essential role, since it yields a dead zone in the optimal 
PRC while a sinusoidal signal does not (see Fig. E]). In other words, oscillators that are 
entrained by stimuli other than insolation may not exhibit a dead zone in their PRCs. This 
is indeed found in mammals. Mammals possess a master clock in their suprachiasmatic 
nucleus (SCN), which receives light stimuli via retinal photoreceptors, and peripheral clocks 



in body cells 



The peripheral oscillators are entrained by several stimuli such as feeding 



and signals from the SCN through chemical pathways (e.g., hormones), and glucocorticoid is 



the suspected synchronizer 



36 



37| . From a study of injections of glucocorticoid, Balsalobre 



et al. [38| reported that the PRCs of the peripheral oscillators in the liver did not have a 
dead zone. 

Our result also agrees with other experimental observations. Our theory implies that a 
dead zone should be located where the concentration £LC,fc(0) * s l° w (0 < < 7r in Fig. EB), 
and that to achieve optimality, the concentration of £lc,/c(0) should be maximal in the region 
where the PRCs exhibit a large phase shift. In Drosophila, the timeless gene is regarded as 
the molecular implementation of £lc,/c(0)- It is experimentally known that light enhances 
the degradation of the gene product (the TIM protein) {^J, and the TIM protein peaks 
during the late evening. Figure [7] shows observations of the PRC of Drosophila against light 
pulses as a function time (hour) from Ref . |40j] ; circles describe the experimental data and the 
solid line expresses a simple moving average (three-points average), respectively. Because 
the center of the part of the PRC that can be phase shifted approximately corresponds 
to the peak of the concentration, as denoted above, when estimated from the PRC alone, 
the concentration peak of the TIM protein should occur at about 18 h. This time is also 
close to the experimental evidence (i.e. late evening). Therefore, our theory can be used to 
hypothesize further molecular behavior affected by the zeitgebers. 

In summary, we have constructed a model that regards circadian oscillators as a global 
optimization of entrainability and regularity. Our model is different from conventional sys- 
tems' biology in the sense that we rely on neither exhaustive computer simulations nor 
species-specific reaction models. We have shown that our model is consistent with much 
experimental evidence as mentioned above. The extension and improvement of our method 
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are 



possible and they are left as an area of future study. 
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FIG. 1: Illustration of typical PRCs (A) with and (B) without a dead zone. Inside a dead zone, 
the input signal neither advances nor delays the clock. 
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FIG. 2: (A) Insolation input of Eq. fT2| where the shaded region denotes night. (B) Time course 
of xlc,/c (0) (Eq. [15]) , which is a variable to be multiplied by the parameter p (Eq. [Fi|) . 
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FIG. 3: (A)-(C) Landscape of *$>(A,5) as a function of A and 5 with insolation input (Eq. I12|) 
for (A) a = 0, (B) a = 0.5, and (C) a = 1.0, where the maximum points are parameters for the 
optimal PRC. (D)-(F) Optimal PRCs with insolation input: (D) a = 0, (E) a = 0.5, and (F) 
a = 1.0. In (D)-(F), the solid and dashed lines denote iPRCs Uj(<fi) and pPRCs Z(<f>), respectively 
(in (D), solid and dashed lines are indistinguishable). The maximal parameters (A, 5) for (D)-(F) 
are (D) (vr,0), (E) (2.31,1.99), and (F) (2.30,2.72). In (D), a parallel shift of the PRC is also 
optimal (5 can be an arbitrary value). In (E), symmetric PRCs with respect to the horizontal axis 
are also optimal. In (F), symmetric PRCs with respect to the horizontal axis or = 37r/2 are also 
optimal (see the text). The pPRCs correspond to experimentally observed PRCs. 
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(A) (B) 




FIG. 4: (A) Dual-axis plot for the dependence on a of the dead-zone length L (solid line; left axis), 
and the entrainability ratios Rjj (dashed line; right axis) and Rz (dot-dashed line; right axis). 
Rjj and Rz are the ratios of the entrainability of the optimal PRC to that of the sinusoidal iPRC 
(Eq.[T7|) and the pPRC (Eq. [T8|) . respectively. (B) Comparison of several pPRCs Z{<j>): the optimal 
pPRC (solid line), the pPRC of the sinusoidal iPRC of Eq. [T7] (dashed line), and the sinusoidal 
pPRC of Eq. [18] (dot-dashed line) for a = 0.5. These results have the same period variance a\ = 1. 




FIG. 5: Dual-axis plot of iPRC Uj((p) (solid line; left axis), pPRC Z{(p) (dashed line; left axis), and 
the concentration ielc fc(0) (dot-dashed line; right axis), for the (A) optimal and (B) alternative 
cases. In (A), the dead zone is located where the concentration of xlc.A;^) is minimal, whereas in 
(B), it is located where the concentration is maximal. Although the iPRCs are the same in both 
cases, the pPRC of the optimal case yields better entrainability. 
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FIG. 6: Optimal PRCs with sinusoidal input (Eq. [T3|) : (A) a = 0.5 and (B) a = 1.0, where the 
solid and dashed lines denote iPRCs Uj(4>) and pPRCs Z((p). For a = 0, the optimal PRCs are 
sinusoidal and are equivalent to those of Fig. [3p. The maximal parameter (A, 5) is (ir, 0) in all 
cases. PRCs that are symmetric with respect to the horizontal axis are also optimal. The pPRCs 
correspond to experimentally observed PRCs. 
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FIG. 7: Experimentally observed PRC in Drosophila with light pulses (circles) and its moving 
average (solid line). Shaded and nonshaded regions indicate subjective night and day, respectively. 



The data were originally obtained by Konopka and Orr, and we took values from Fig. 2 in Ref. 
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Supporting information: 
Circadian clocks optimally adapt to sunlight for reliable 

synchronization 

Yoshihiko Hasegawa and Masanori Arita 



This supporting information describes the detailed calculations for obtaining the optimal PRCs for 
circadian oscillators, including procedures for deriving regularity and entrainability measures. We also 
present optimal PRCs derived from another insolation input (the squared sinusoidal insolation signal) . 
Equation and figure numbers in this section are prefixed with S (e.g., Eq. SI or Fig. SI). Numbers 
without the prefix (e.g., Eq. 1 or Fig. 1) refer to those in the main text. 

1 Definition of regularity 

We show a detailed procedure for calculating regularity. In the presence of molecular noise, the period 
varies cycle to cycle as illustrated in Fig. IS1IA. where Tj denotes a period of i-th cycle. We quantify 
regularity by the variance of the period variation due to the stochasticity. 

We can naturally define phase <j> <G [0, 2ir) on the unperturbed oscillation of Eq. 1 by 

I- - < S1 > 

where f2 is the angular frequency of the oscillation. The phase <j> in Eq. [Si] is only defined on a closed 
orbit of the unperturbed limit-cycle oscillation. However, we can expand the definition into the entire 
x = (xi,X2,--- ,xjv) space, where the equiphase surface is referred to as the isochron I(<j>) (Fig. 15TB). 
A phase-reduction approach pQ reduces the perturbed iV-dimensional dynamics into a one-dimensional 
differential equation with respect to the phase <f> defined above. By using stochastic phase reduction, 
the iV-dimensional limit-cycle oscillator subject to noise (Eq. 2) can be transformed into Langevin 
equation 3: 

i, n 

^=fi + X;W)Q<(06(i), (S2) 

i=l 

where Ui{4>) is the iPRC, Qi{4>) is a multiplicative noise term and is white Gaussian noise with 
the correlation — 28ij8(t — t 1 ) (see the main text). Let P(<f>;t) be the probability density 

function of <f> at time t. From Eq. IS2| the Fokker-Planck equation (FPE) of P ((/)', t) is given by 

where 

N , 

N 

g{4>) = X>(0 a QiM a . (S5) 



i 



By a slow variable ip — <f)~ fit introduced in the main text, the FPE of the probability density function 
IL(tp;t) = P(4> = ip + Qt;t) is 

Q- t n&;t) = \-Q^Hf + nt) + g^g& + m) j i% ; t). (se) 

With sufficiently weak noise, H(<p; t) is a slowly fluctuating function of t. In such cases, F(ip + fit) 
and Q(ip + fit) fluctuate much faster than H(ip; t), thus these two terms can be averaged for one period 
while keeping H((p; t) constant (phase averaging). By phase averaging, J-(ip + fit) vanishes because of 
the periodicity (use integration by parts): 

§- t n(v;t) = D-^ii(<p-,t), (S7) 



i r 2rr N 

D = ^ deY / U l (9) 2 Q t (9) 2 . (S8) 



with 



i=l 

Please see Ref. [1 for further details of stochastic phase reduction and the phase-averaging procedure. 
From Eq. IS71 P(4>; t) = H(ip = <fi — fit; t) obeys a simple one-dimensional diffusion equation, and its 
solution is represented by 

which shows that the variance of the phase after one period is = 2DT. We approximate the period 
variance Vt, which is the variance of the first passage time from to 27r, with the phase variance V^, 
after Kori et. al [2]- As the phase cj> crosses a threshold <fi — 2tt with gradient 2ir/T without noise, 
there is a scaling relation *JVt — \/^4>T/(2tt) for sufficiently weak noise (Fig. IS1C). Consequently, the 
variance of the period is approximated by 

v T ~ 2dt (L^ 2 = J 2 J d e J2 u^efQrio) 2 , (sio) 

which is Eq. 4. 

In order to see the reliability of Eq. ISIOl we carried out Monte Carlo (MC) simulations of a 
stochastic limit-cycle oscillator (we followed [3] for the MC algorithm). We employed a typical genetic 
oscillator model (a smooth oscillator) from Novak and Tyson [J] ; it is given by the following Langevin 
equation with respect to the molecular species concentrations [x\, X2, X3): 



~~dT = K p + x p ~ + VQW'' 



*E1 - k x,-u, X2 

j, — r^sy^x v ay 



dt ° y y K m + x 2 J 

dx 3 

—7- = K SZ X 2 - KdzX 3 , 

at 

where V sx = 5.44, k dx = 8.99, k sz = 4.49, k dz = 4.49, k sy = 18.0, k dy = 0.0, V dy = 2.72, K m = 0.00303, 
K d = 0.303, and p — 4.0 (we employed the same parameter values as in Ref. [5]; the oscillator has a 
unit period T = 1 with these parameters). £(i) is white Gaussian noise with correlation (£{t)£(t')) = 
26 (t — t'); q represents the noise intensity (also see the main text). The period is calculated as the 
difference in the time between the first and the second passages when x 2 crosses 0.5 from below. We 
discarded those with periods smaller than 0.2 to remove very small values that were due to a stochastic 
threshold crossing. Figure [Sit) shows the standard deviation of the period (i.e., \/Vt) as a function 
of the noise intensity q, where the solid line and circles represent the values of Eq. IS 101 and those of 
the MC simulation, respectively. We see excellent agreement between the analytic and MC results for 
q < 10 ~ 3 , which verifies the reliability of Eq. IS 101 
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Figure SI: (A) Time course of a hypothetical circadian oscillator subject to molecular noise, where the 
period Ti varies cycle to cycle. (B) Illustration of the isochron. The solid and dashed lines describe 
a limit-cycle trajectory and its isochron drawn at intervals of 7r/6, respectively. (C) Relation between 
the phase variance V^ and the period variance Vt in Langevin equation IS2I (the solid lines represent 
trajectories of the Langevin equation). is the variance of the phase <fi a t time t = T and Vt is the 
variance of the first passage time from to 2tt, which can be approximated by Vt — V<j>T 2 / (2tt) 2 . (D) 
Standard deviation of the period y/Vr as a function of the noise intensity q, in a log-log scale. The 
solid line describes the analytic results from Eg. ISlOl and the circles show those of the MC simulations. 
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Figure S2: (A) Illustration of the existence of a stable solution in ip = Att + X©(^) = 0. Solid and 
dashed lines represent x0(V>) and — Af2, respectively, and an intersection point with d®/dip < is a 
stable solution (an open circle). (B) Arnold tongue (colored region), which shows the parameter region 
for synchronization to an external signal, with respect to the signal angular frequency oj (vertical axis) 
and the signal strength x (horizontal axis) . The dashed line is a linear approximation of the border of 
the Arnold tongue when the input strength \ is sufficiently small. 



2 Definition of entrainability 

We show a detailed procedure for calculating entrainability. Considering the parameter perturbation 
p + dp due to the light stimuli, we obtain Eq. 2 

= Fi(x;p + dp), 

where we do not take into account the effects of noise, as explained in the main text (Qi{x) — 0). 
Because we assume that the perturbative term dp is sufficiently small, we can approximate the dynamics 
of the phase <f> £ [0, 2ir) with the phase-reduction approach, as follows: 

dt 4^ x dxi dp 

= n + x z(<j>) P (iot), (sii) 

where Z(<f>) is a pPRC (Eq. 5). The dynamic equation for the slow variable ip — <fi — uit is 

^ = AO + X e(V), (S12) 

where Afl = O — u> and <d(ijj) are defined in Eq. 6. When ip = in Eq. IS 121 has a stable stationary 
solution, the oscillator can synchronize to a periodic input signal p(u>t). A stable solution of if) = 
is an intersection point of and —Aft with dQ/dip < (an empty circle in Fig. IS2A). Then a 

condition for the existence of a stable solution is 

xe(v m ) + n<w< x e(^) + a (si3) 

where tpM = argmax^0(-0) and tp m — argmin^,0(^). The extent of synchronization is quantified by 
a domain known as the Arnold tongue, which shows the synchronization condition with respect to \ 
(the signal strength) and u) (the angular frequency of the input signal) . The shaded region in Fig. IS2B 
represents the Arnold tongue as functions of \ an d w; inside the tongue the oscillator of interest can 
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synchronize to a periodic input signal. Equation IS 131 constitutes a linear approximation of the Arnold 
tongue for sufficiently small x- The extent of synchronization, the entrainability, can be measured by 
the width of the Arnold tongue, which is given by \ (©(V'm) — ©(VvO) under the linear approximation. 
In the present paper, we define the entrainability £ as a proportionality coefficient of the width (cf. 
Eq. 7). 



3 Variational method 

From Eq. 8, the variational equation to be optimized reads 



C[U] 



£[U] - XV T [U], 



1 

2Wo 



2tt N 



W-M-Pie-MWe) 



dFi{6;p) XT 3 



dp 2tt 2 

where A is the Lagrange multiplier. The variational condition 5C[U] = yields iPRC 

7T 2 p((f> - 1p M ) - P{4> - 0m) dFi{4>\ p) 



Ui{6fQi{ef 



Ui{<j>) 



T 3 A 



>W>) 2 



dp 



and pPRC 



jf_ IpM ) P(<f> ~1pm) ( 9F% (0; P) 

T 3 A 



Qi(4>? 



V dp 



Letting the variance of the period be Vt — Oj~, the Lagrange multiplier is 



\ 4T 3 4 J 



V dp 



Substituting Eas. I5l6l and I3T71 into Eq. 7, we have £(A,8) in Eq. 10. 



,(S14) 



(S15) 



(S16) 



(S17) 



4 Explicit expression of ^(A,^) 

Insolation input 

For the insolation input of Eq. 12, ^(A, 6) (Eq. 11) is given by 



*(A,J) 



* a (A,(5) 0<A<7T, 
*b(A,5) 7T<A<27T, 



with 



(S18) 



* a (A,<5) = [-3a 2 sin(3A + 28) + 32a cos(2A + 8) + 12a 2 A cos(2<5 + A) 

-6a 2 sin(2<5 + A) + 12a 2 7r cos(A + S) 2 - 128a cos(A + 6) 

+ {-24a 2 7rcos((5) 2 + (128a + 18a 2 sin 8) cos 8 + (-12tt + 24A)a 2 - 48tt + 48A} cos A 



+12 ( - sin A + 7T J a 2 cos((5) 2 + (24a 2 7r sin A sin 8 - 32a) cos 8 
+ (-64a sin 8 - 48 - 27a 2 ) sin A + 48?r + 12a 2 7r] , 



where \&b(A,5) = 4 r a (— A + 27r, — 5 + 2?r). Equation IS 181 and the corresponding optimal PRCs were 
plotted in Fig. 3 for a = 0, 0.5, and 1.0. 
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(A)a = 




Figure S3: (A)-(C) Landscape of '■['(A,^ (cf. Eg. IS19[) as functions of A and 6 with sinusoidal input 
(Eq. 13) for (A) a = 0, (B) a = 0.5, and (C) a — 1.0, where the maximum points are parameters for 
the optimal PRC. 



Sinusoidal input 

For the sinusoidal input of Eq. 13, \I/(A,<5) is given by 



*(A, S) = — (1 - cos A) [-a 2 cos(26 + A) + 2a 2 + 4l 
2q 



(S19) 



We plot Eq. lS19l as functions of A and 8 in Fig. IS3[ where the optimal parameter A is it regardless of 
a. The corresponding PRCs for a — 0.5 and 1.0 have been plotted in Fig. 6 (for a = 0, the optimal 
PRCs are sinusoidal). 



5 Sinusoidal iPRC and pPRC 

Sinusoidal iPRC 

An explicit expression of Eq. 17 (sinusoidal iPRC) is 



/47r 2 CT 2 

U 3 (fa C ) = \J Shl + C ) ' 

which yields the period variance of Vt = of. Then the corresponding pPRC is given by 



Z(6:c) 



Uw 2 al . ,, . ,„ 
— sin (0 + c) (1 — a sin < 



(S20) 



(S21) 



where we used Eq. 5. 
Sinusoidal pPRC 

For the pPRC Z(<p) to be a sinusoidal function, the iPRC Uj{<j>) must be 



dF 2 {^_p)Y 1 
dp 



sin(0 + c), 



where we used Eq. 5. An explicit expression of Eq. IS22I is 

1 sin(0 + c) 



Uj(fac) = 



Af(c) 1 — a sh 



(S22) 



(S23) 
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Figure S4: Squared sinusoidal insolation input for Eg. IS271 where the shaded region denotes night, 
where Af(c) is a normalizing term 
Af(c) = 



qT 3 a 2 -{(2vT — a 2 — 3) a 2 — 2V1 - a 2 + 2} cos(2c) 



a 2( 1 _ a 2)3/2 



Equation IS23I is normalized so that the period variance becomes Vt = erf,. Using Eq. 5, the corre- 
sponding pPRC is a sinusoidal function: 



(S24) 



which is an explicit expression of Eq. 18. 



Entrainability calculation 



For both Eqs. IS21I and IS241 it is necessary to determine a parameter c that yields the maximum 
entrainability. From Eq. 7, the entrainability for Eqs. IS2ll and IS241 is 



S(A, §,c) = ± [ d6 Z(6 + 5; c) {p(9 - A) - p(9)} , 
27r Jo 

where p(9) is the insolation input. The entrainability is then calculated by 



£ = max£(A, 5, c), 

A,5,c 



(S25) 



(S26) 



which was optimized by the differential evolution 6 . The optimal PRCs are obtained by substituting 
ipM = A + 5 and if) m = 5 into Eqs. IS15I and IS161 



6 Another insolation input signal 

In the main text, we considered a truncated sinusoidal function (Eq. 12) for the insolation input. In 
order to see the effect of the shape of the insolation input, we here use a squared sinusoidal function: 



p(ojt) = ramp(sin(wt)) 2 , 



(S27) 



where the ramp function being defined as ramp(x) = x for x > and ramp(i) = for x < 0. The 
actual insolation does not strictly obey the cosine law due to the effect of diffusion in the atmosphere, 
and the insolation smoothly varies as a function of time. The squared sinusoidal insolation input 
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Figure S5: (A)-(C) Landscape of ^(AjJ) as a function of A and 8 with squared sinusoidal insolation 
input (Fig. IS4I and Eq. IS27|) for (A) a = 0, (B) a = 0.5, and (C) a — 1.0, where the maximum points 
are parameters for the optimal PRC. (D)-(F) Optimal PRCs with squared sinusoidal insolation input: 
(D) a = 0, (E) a = 0.5, and (F) a = 1.0. In (D)-(F), the solid and dashed lines denote iPRCs Uj(<j>) 
and pPRCs Z(<j>). The maximal parameters (A, 8) for (D)-(F) are (D) (tt,0), (E) (1.93,2.18), and (F) 
(1.81,2.24). In (D), a parallel shift of the PRC is also optimal (8 can be an arbitrary value). In (E) 
and (F), symmetric PRCs with respect to the horizontal axis are also optimal. 

of Eq. IS27I is smooth at the boundaries of and ir, and thus can be considered as the insolation 
approximation reflecting the diffusion effect. Figure [S4l shows the squared sinusoidal insolation input 
of Eq. IS27| where the shaded and nonshaded regions represent night and day, respectively. 

Figures [S5K-C portray the landscape of ^(A, 8) as functions of A and 8, and Figs. IS5D-F portray 
the optimal iPRCs Uj{4>) and pPRCs Z{<f) for the squared sinusoidal insolation input (we employed 
T = 1, a% = 1, and q = 1 as in the main text). Figures IS5K-C show \&(A, 8) for (A) a = 0, (B) 
a = 0.5, and (C) a = 1.0, and Figs. IS5D-F show the optimal iPRCs (the solid line) and pPRCs (the 
dashed line) for (D) a = 0, (E) a = 0.5, and (F) a = 1.0. As in the normal insolation input, the 
optimal PRC of a = does not have a dead zone, and its shape is close to the sinusoidal function. 
For a = 0.5, the optimal parameters are (A, 8) — (1.93,2.18) and (4.36,4.10), which yield the dead 
zone in the optimal PRC (Fig. IS5E) because A ^ it. From Fig. IS5F, we see that the optimal PRC of 
a = 1.0 also has a dead zone (parameters are (A, 8) — (1.81,2.24) and (4.47,4.05)). For the extent 
to which the input signal affects the kinetic rates in chemical reactions, we observe qualitatively the 
same optimal PRCs for both types of insolation signals. 
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